#### Cosponsorship Community Detection #####
source("/Users/brucedesmarais/Dropbox/SenPE/Code/CreateSenPENets.R")

## Make the lists of cosponsorship matrices, sponsorship counts, and member ids

senmats <- list()
senspons <- list()
senmems <- list()
for(i in 93:108){
	if(i < 100){
		matfile <- paste("/users/brucedesmarais/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_matrices/senate_matrices/0",i,"_senmatrix.txt",sep="")
		idfile <- paste("/users/brucedesmarais/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_members/senate_members/0",i,"_senators.txt",sep="") 
	}
	if(i >= 100){
		matfile <- paste("/users/brucedesmarais/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_matrices/senate_matrices/",i,"_senmatrix.txt",sep="")
		idfile <- paste("/users/brucedesmarais/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/senate_members/senate_members/",i,"_senators.txt",sep="")
	}
	senmati <- read.csv(matfile,header=F)
	senmati <- (senmati==1) + (senmati == 2)
	senmems[[i-92]] <- read.csv(idfile,header=F)
	amat <- matrix(0,nrow(senmati),nrow(senmati))
	sponn <- numeric(nrow(senmati))
	for(j in 1:nrow(amat)){
		for(k in (1:nrow(amat))[-j]){
			amat[k,j] <- sum(senmati[j,]*senmati[k,])
			amat[j,k] <- sum(senmati[j,]*senmati[k,])
			
			}
		}
		print(i)
	senmats[[i-92]] <- amat
	senspons[[i-92]] <- sponn
}


# Make list of node-level data

sendat <- read.csv("/users/brucedesmarais/Desktop/Research/TERGM/Data/Cosponsorship_Data_2009_Update/Cosponsorship_Data_2009_Update/SH.csv",header=T,stringsAsFactors=F)
sendats <- list()
for(i in 93:108){
	sendati <- subset(sendat,is.element(sendat$congress,i))
	sendats[[i-92]] <- sendati[match(senmems[[i-92]][,3],sendati$ids),]
}

### Community Detection ###

cosp_comm <- list()
library(igraph)

for(t in 4:13){
g2 <- graph.adjacency(senmats[[t]],weighted=T,mode="undirected")
cg2w <- spinglass.community(g2)
memw <- cg2w$membership
cosp_comm[[t]] <- list(cg2w,memw)
}

com_mod <- NULL

fowler <- read.csv("SH.csv",stringsAsFactors=F)
fowler <- subset(fowler,fowler[,2]=="S")


for(t in 4:13){
g <- graph.adjacency(senmats[[t]],weighted=T,mode="undirected")
mem <- cosp_comm[[t]][[1]]$membership
maxMod <- modularity(g,mem,weight=E(g)$weight)
senmemi <- as.character(senmems[[t]][,1])
congi <- congs[t]
idmi <- subset(id_master,id_master$Congress==congi)
sum_idi <- idmi$cosp_summary[match(sub("  "," ",senmemi),idmi$cosp_senmem)]
fowi <- subset(fowler,fowler$congress==congi)
party <- sendats[[t]]$party
party <- as.numeric(party==200)
party[which(is.na(party))] <- 0
partyMod <- modularity(g,party,weight=E(g)$weight)
com_mod <- rbind(com_mod,c(congi,maxMod,partyMod))

}

# plot modularity vs. polarization in cosponsorship
pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/com_mod.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(com_mod[,c(1,2)],ylim=c(min(c(com_mod[,c(2,3)])),max(c(com_mod[,c(2,3)]))), lwd=2.5,xlab="",ylab="Modularity",type="l",yaxt="n")
lines(com_mod[,c(1,3)],col="grey55",lwd=2.5)
abline(v=96:105,lwd=2.5,col="grey70",lty=3)
points(com_mod[,c(1,2)],cex=1.5,pch=19)
points(com_mod[,c(1,3)],pch=23,cex=1.5,col="grey55",bg="grey60")
legend("bottomright",legend=c("Max Modularity","Polarization"),col=c("black","grey55"),lwd=2.5,pch=c(19,23),pt.bg=c("black","grey55"),bg="white")
axis(1,at=96:105)
axis(2,at=(3:10)/100)
dev.off()

# Plot colored community membership for cosponsorship

library(colorRamps)
pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/cosp_comm.pdf",height=3.5, width=6, family="Times", pointsize=13.5)
par(las=1,mar=c(3.5,5,.5,.5),cex.lab=1.25)
plot(96:105,seq(0,100,length=length(4:13)),col="white",xlab="",ylab="Cumulative Percent",bty="n",xlim=c(95.5,105.5),xaxt="n")
red2blues <- gray((100:1)/100)

for(i in 4:13){
	memsi <- colnames(mats[[i]])
	csizei <- cosp_comm[[i]][[1]]$csize
	cmemi <- cosp_comm[[i]][[1]]$membership
	coms <- which(csizei>1)-1
	congi <- congs[i]
	csizei <- csizei[which(csizei>1)]
	props <- 100*csizei/sum(csizei)
	pdem <- numeric(length(coms))
	for(j in 1:length(pdem)){
		cmemj <- memsi[which(cmemi==coms[j])]
		party <- sendats[[i]]$party[which(cmemi==coms[j])]
		pdem[j] <- mean(party==100,na.rm=T)
		}
	cszsrt <- csizei[order(pdem)]
	pdem <- pdem[order(pdem)]
	cols <- red2blues[1+round(pdem*99)]
	cssz <- cumsum(cszsrt/sum(cszsrt))
	lb <- c(0,cssz[1:(length(cssz)-1)])*100
	ub <- cssz*100
	for(j in 1:length(pdem)){
		rect(xleft=congi-.5,xright=congi+.5,ybottom=lb[j],ytop=ub[j],col=cols[j])
		}
	}

axis(1,at=96:105)
dev.off()

